home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / ASTRNOMY / AA_51.ZIP / REFRAC.C < prev    next >
C/C++ Source or Header  |  1987-12-02  |  2KB  |  69 lines
  1. /* Atmospheric refraction
  2.  * Returns correction in degrees to be added to true altitude
  3.  * to obtain apparent altitude.
  4.  *
  5.  * -- S. L. Moshier
  6.  */
  7.  
  8. extern double atpress; /* millibars */
  9. extern double attemp; /* degrees C */
  10. extern double DTR; /* pi/180 */
  11.  
  12. double refrac(alt)
  13. double alt;    /* altitude in degrees */
  14. {
  15. double y, y0, D0, N, D, P, Q;
  16. int i;
  17. double tan();
  18.  
  19. if( (alt < -2.0) || (alt >= 90.0) )
  20.     return(0.0);
  21.  
  22. /* For high altitude angle, AA page B61
  23.  * Accuracy "usually about 0.1' ".
  24.  */
  25. if( alt > 15.0 )
  26.     {
  27.     D = 0.00452*atpress/((273.0+attemp)*tan( DTR*alt ));
  28.     return(D);
  29.     }
  30.  
  31. /* Formula for low altitude is from the Almanac for Computers.
  32.  * It gives the correction for observed altitude, so has
  33.  * to be inverted numerically to get the observed from the true.
  34.  * Accuracy about 0.2' for -20C < T < +40C and 970mb < P < 1050mb.
  35.  */
  36.  
  37. /* Start iteration assuming correction = 0
  38.  */
  39. y = alt;
  40. D = 0.0;
  41. /* Invert Almanac for Computers formula numerically
  42.  */
  43. P = (atpress - 80.0)/930.0;
  44. Q = 4.8e-3 * (attemp - 10.0);
  45. y0 = y;
  46. D0 = D;
  47.  
  48. for( i=0; i<4; i++ )
  49.     {
  50.     N = y + (7.31/(y+4.4));
  51.     N = 1.0/tan(DTR*N);
  52.     D = N*P/(60.0 + Q * (N + 39.0));
  53.     N = y - y0;
  54.     y0 = D - D0 - N; /* denominator of derivative */
  55.  
  56.     if( (N != 0.0) && (y0 != 0.0) )
  57. /* Newton iteration with numerically estimated derivative */
  58.         N = y - N*(alt + D - y)/y0;
  59.     else
  60. /* Can't do it on first pass */
  61.         N = alt + D;
  62.  
  63.     y0 = y;
  64.     D0 = D;
  65.     y = N;
  66.     }
  67. return( D );
  68. }
  69.